Genetic–environment associations explain genetic differentiation and variation between western and eastern North Pacific rhinoceros auklet (Cerorhinca monocerata) breeding colonies

Abstract Animals are strongly connected to the environments they live in and may become adapted to local environments. Examining genetic–environment associations of key indicator species, like seabirds, provides greater insights into the forces that drive evolution in marine systems. Here we examined a RADseq dataset of 19,213 SNPs for 99 rhinoceros auklets (Cerorhinca monocerata) from five western Pacific and 10 eastern Pacific breeding colonies. We used partial redundancy analyses to identify candidate adaptive loci and to quantify the effects of environmental variation on population genetic structure. We identified 262 candidate adaptive loci, which accounted for 3.0% of the observed genetic variation among western Pacific and eastern Pacific breeding colonies. Genetic variation was more strongly associated with pH and maximum current velocity, than maximum sea surface temperature. Genetic–environment associations explain genetic differences between western and eastern Pacific populations; however, genetic variation within the western and eastern Pacific Ocean populations appears to follow a pattern of isolation‐by‐distance. This study represents a first to quantify the relationship between environmental and genetic variation for this widely distributed marine species and provides greater insights into the evolutionary forces that act on marine species.


| INTRODUC TI ON
Examining questions about local adaptation for vagile marine species that have broad distributions is of interest because populations living in different oceanographic regions are exposed to different ecological and environmental conditions.Therefore, ecological adaptation to these environments can lead to reproductive isolation and genetically distinct populations (Friesen, 2015;Friesen et al., 2007;Munro & Burg, 2017).Using genomic analyses, we can build on studies of marine species that used traditional markers (microsatellites and mitochondrial DNA; reviewed in Taylor & Friesen, 2012) to delineate population boundaries and characterize population structure and test for associations between genetic and environmental variation.
Examining the relationship between environment and genetic variation for marine species will provide greater insights into the forces generating diversification.Determining the role that environmental variation plays on genetic variation is important (Garcia De Leaniz et al., 2007;Miller et al., 2018;Razgour et al., 2019) given the rapid and substantial climate changes occurring in marine environments (Carozza et al., 2019;Constable et al., 2014;Doney et al., 2012).
Studying genetic-environment associations will help to predict how adaptable populations are and how they will respond to changing climates (Capblancq et al., 2020), given that changes in environment can affect fitness (Jenouvrier et al., 2018;Sydeman et al., 2012).
Seabirds offer a compelling system to examine questions about adaptative genetic variation in marine environments because their population dynamics are linked with environmental change and variation (Piatt et al., 2007).For these reasons, seabirds are considered key indicator species of marine systems.Seabirds often exhibit strong natal and breeding philopatry, but philopatry can be quite variable between species (Coulson, 2016).If adult and juvenile philopatry are high, there is potential for local adaptation at relatively small spatial scales (Stiebens et al., 2013) and population genetic patterns may reflect local adaptation.If philopatry is low, however, then the potential for individuals to disperse from their breeding and natal colonies may increase gene flow among colonies (Castillo-Guerrero et al., 2020;Hipfner et al., 2020;Quillfeldt et al., 2017).As a result, individuals will be exposed to a wider range of environmental conditions and thereby reduce the potential for local adaptation (Antoniou et al., 2021;Jahnke et al., 2022;Lombal et al., 2020).This high dispersal potential is reflected in the complex population genetic patterns found for many species (Burg & Croxall, 2001, 2004;Dhami et al., 2013;González-Jaramillo & Rocha-Olivares, 2011;Hall et al., 2009;Milot et al., 2008;Munro & Burg, 2017).In this study, we examine the relationship between genetic and environment variation in the rhinoceros auklet (Cerorhinca monocerata), a colonial nesting seabird broadly distributed across the northern Pacific Ocean.This species forages on zooplankton and small schooling fish, which rhinoceros auklets capture underwater during shallow dives (mean dive depth ranges between 9.2 ± 6.7 m and 12.1 ± 5.5 m; Cunningham et al., 2018;Kato et al., 2003).Rhinoceros auklets use both offshore and coastal waters to forage for food depending on the dominant foraging fish community (Iida et al., 2024).
F I G U R E 1 (a) Range map of the rhinoceros auklet (Cerorhinca monocerata) across the northern Pacific Ocean.Yellow represents resident areas, while blue represents wintering areas.Shapes represent the 15 breeding colonies where samples were collected from (b) the five breeding colonies sampled from Japan; (c) the three breeding colonies sampled from coastal Alaska; (d) the four colonies sampled from coastal British Columbia; and (e) the three breeding colonies sampled from coastal Washington and coastal California.Breeding colony abbreviations match those found in Table 1.Map was created using QGIS 3.12 with range map layers provided by Bird Life International and Handbokok of the Birds of the World (2022, v 2022.1).

(a) (b) (c) (d) (e)
Previous studies using microsatellites have shown that populations on either side of the Pacific Ocean are genetically distinct (Abbott et al., 2014;Hipfner et al., 2020).Eastern Pacific Ocean rhinoceros auklets are also smaller than rhinoceros auklets from the western Pacific Ocean (Hipfner et al., 2020).Western and eastern populations forage for different forage fish (Burger et al., 1993;Cunningham et al., 2018;Takahashi et al., 2001;Wilson & Manuwai, 1986) and thereby exhibit foraging behavior and diet differences (Cunningham et al., 2018;Senzaki et al., 2014;Takahashi et al., 2001).Finally, the wintering grounds of eastern and western populations (Figure 1) are also distinct (Hipfner et al., 2020), which may further contribute to the genetic differences previously reported.As oceanic conditions continue to change, in response to climate change, and extreme oceanic events like marine heatwaves (Oliver et al., 2018) become more prominent, questions remain about how climate change will affect marine species like rhinoceros auklets (Wagner et al., 2023).Given that rhinoceros auklets have a broad distribution across the North Pacific, examining how genetic variation is associated with environmental variation in this species can provide greater insight into whether populations are adapted to their conditions and how increasing climate variability will affect populations across their distribution.We used a restriction-site associated DNA (RADseq) dataset to examine whether genetic differences between western Pacific and eastern Pacific populations are associated with oceanic environmental variation at breeding colonies.
Additionally, we examined genetic variation within western Pacific and eastern Pacific populations to determine whether breeding colonies are genetically distinct and whether adaptive genetic variation influences population genetic differences.This last question is especially relevant because rhinoceros auklets breed across a wide geographic area in the northern Pacific Ocean (Figure 1).We hypothesize that sea surface temperature will account for the largest component of genetic variation out of the three environment variables because temperature varies between the eastern and western Pacific Ocean and along a latitudinal gradient within each basin (whereby sea surface temperatures are warmer at southern breeding colonies and cooler at more northern temperatures).Sea surface temperature is a strong predictor of survival and reproduction rates (Jenouvrier et al., 2018;Sydeman et al., 2012) and therefore has a strong effect on seabird fitness.Given that genetic variation and structure in marine systems is associated with population connectivity (i.e., spatial separation during the breeding season; Cowen et al., 2007;Friesen et al., 2007), we evaluate the roles of geographic distance and population connectivity (using least-cost path data) on genetic variation.Examining multiple hypotheses about the source of genetic variation observed for this widely distributed species will provide greater context and insights into how climate and environmental variation affect genetic variation for marine species.

| DNA extraction and library preparation
DNA was extracted from blood samples using a salting out extraction protocol (for samples from the eastern Pacific, Miller et al., 1988) or a Qiagen DNeasy kit (for samples from western Pacific).Genomic DNA was used to construct nextRAD genotyping-by-sequencing libraries (SNPsaurus, LLC) using the Sbf1 enzyme as described by Baird et al. (2008).Genomic DNA was first fragmented with Nextera reagent (Illumina, Inc), which also ligates short adapter sequences to the ends of the fragments.The Nextera reaction was scaled for fragmenting 15 ng of genomic DNA, although 20 ng of genomic DNA was used for input to compensate for the amount of degraded DNA in the samples and to increase fragment sizes.Fragmented DNA was then amplified for 27 cycles at an annealing temperature of 74°C, with one of the primers matching the adapter and extending 10 nucleotides into the genomic DNA with the selective sequence GTGTAGAGCC.Only those fragments starting with that sequence can be hybridized by the selective sequence of the primer and efficiently amplified.This protocol resulted in a final library fragment size of 450 bp (Etter et al., 2011).The nextRAD libraries were sequenced on an Illumina NovaSeq 6000 with one lane of single-end 150 bp reads.All genomic library preparations and sequencing were completed at the University of Oregon.
Sequences were demultiplexed and then trimmed to 122 bp by SNPsaurus using the SNPsaurus pipeline with the bbduk package (BBMap tools, http:// sourc eforge.net/ proje cts/ bbmap/ ).Next, we assembled reference loci by collecting 10 million high-quality reads, evenly from all of the samples (~70,000 reads per individual were used), and excluding loci with fewer than seven or more than 700 reads.This range of seven to 700 represents a standardized number calculated by SNPsaurus to retain as many loci as possible without compromising the quality of the data with low-quality reads.Overall mean depth of the reference genome was 65×.Loci that met the previously stated criteria were then aligned to the assembled reference genome using custom script from SNPsaurus (SNPsaurus, LLC).
For the de novo alignment, we mapped 152,204,819 of the original 289,864,865 single-end reads to the de novo reference genome using an identity threshold of 0.95 using bbmap (BBMap tools).
Genotype calling was done using the callvariants tool (BBMap tools), with the following settings (multisample = t rarity = 0.05 minallelefraction = 0.05 usebias = f ow = t nopassdot = f minedistmax = 5 minedist = 5 minavgmapq = 15 minreadmapq = 15 minstrandratio = 0.0 strandedcov = t).The genotype data were converted to a VCF file where we filtered the data to remove loci with a minimum frequency of less than 3%, a Q-score below 20 and removed all individuals with greater than 60% missing data (an additional 13 individuals did not meet this criterion and were omitted from all analyses).The average percentage of missing data was much lower than this original threshold (mean = 5.8% missing data; median = 3.3% missing data), although we included three individuals with 40% missing data because they grouped with other individuals from the same population.To ensure that relatedness did not skew our results, we calculated relatedness among individuals in Genodive 3.04 (Meirmans, 2020).Relatedness among individuals from the same population was <0.08, with exception to one pair that had a relatedness of 0.26, suggesting that one set of full siblings from Teuri were present in our data.Given the low level of relatedness among our data, we retained all samples in our analyses.We retained all of the 19,213 SNPs following the filtering for our examination of genetic-environment associations.

| Population genetic structure
To characterize genetic population structure, we filtered our full SNP dataset to allow only a single SNP per contig because SNPs on the same locus are physically linked.SNPs were also checked for linkage disequilibrium in PLINK (version 1.9; Purcell et al., 2007) to further ensure that we did not include any loci that were correlated with each other.To characterize population genetic structure, we calculated Cavalli-Sforza and Edwards distance (Cavalli-Sforza & Edwards, 1967) between all individuals and ran a principal coordinate analysis (PCoA) on the distance matrix in PAST 4.11 (Hammer et al., 2001).Additionally, we examined population genetic structure with discriminant analysis of principal components (DAPC) using the adegenet package (Jombart, 2008) in R 4.2.2.The purpose of using this analysis was to provide an alternative method to examine population genetic support to our PCoA.DAPC converts the raw genetic data to principal components and then uses discriminant function analysis to detect population structure.DAPC is considered more sensitive for detecting population structure than PCoA because it reduces within group variation while optimizing between group variation.We ran k-means clustering analysis from 1 to 15 (based on the number of sampled breeding colonies) and selected the optimal K by choosing the K with the lowest Bayesian information criterion (BIC; Jombart, 2008).Finally, we calculated pairwise F ST values between all 15 breeding colonies to characterize the level of genetic differentiation across the range in GENODIVE; deviations from zero were determined using 10,000 permutations.
We downloaded 21 ocean environmental data variables from the Bio-Oracle database Version 2.2 (https:// www.bio-oracle.org/ ; et al., 2018).This dataset contains marine data layers for ecologically relevant variables for present day conditions; variables are based on seasonal and monthly averages calculated from the period of 2000-2014.Variables within this dataset have a spatial resolution of 5 arcmin (Assis et al., 2018).Environmental layers were loaded into QGIS 3.28.0(http:// qgis.osgeo.org), and then, we used the point sampling tool to collect environmental data for each breeding colony using the latitude and longitude of each breeding colony as our sampling location.Although rhinoceros auklets breed on islands, the data used represent the marine conditions from the cell that the geographical location falls in.We then tested for intercorrelations between the 21 variables using a Spearman's r correlation.Most variables exhibited correlations that exceeded 0.7, which left us with three variables to use in our analysis of genetic and environment We chose these variables because they are biologically relevant to a variety of marine taxa (including whales, urchins, and fish) and are variables that have been analyzed in previous studies as they are thought to influence habitat occupancy, dispersal, and thereby population genetic structure (Antoniou et al., 2021;Banks et al., 2007;Fullard et al., 2000;Lombal et al., 2020).We used these three variables to characterize environmental variation of ocean conditions at breeding colonies because (a) rhinoceros auklets use coastal waters for foraging (Iida et al., 2024); (b) tracking data for foraging trips outside of coastal waters were not available, and (c) we did not have data available for the wintering areas used by the individuals that we collected genetic samples for.Sea surface temperature is likely an important variable for seabirds, given that a review examining the impact of climate change on seabird populations revealed that sea surface temperature negatively affects survival and reproductive success of seabirds in 90% and 70% respectively of all studies reviewed (Sydeman et al., 2012).Less is known about the effect of pH on seabirds, but changes in salinity can lead to increases in energetic costs (Constable et al., 2014;Gutiérrez, 2014) which may further affect fitness.Ocean acidification alters pH and is predicted to reduce coral populations (Hoegh-Guldberg et al., 2007) which will thereby affect many prey species that depend on the resources provided by these habitats (Chambers et al., 2011;Constable et al., 2014), and has the potential to exhibit large negative effects on the fitness of seabirds, especially if the forage fish they rely on become less abundant and harder to find (Chambers et al., 2011).Finally, oceanic current velocity influences upwelling patterns, which transports nutrients closer to the surface, thereby bringing forage fish closer to the surface.In the absence of upwelling, forage fish remain in deeper waters, forcing individuals to dive deeper for food thereby increasing energetic demands (Soldatini et al., 2023).Although the effect of each environment variable individually affects genetic patterns remains unknown, the combined effect of these factors will affect fitness, a key component of genetic variation and population structure.Additionally, these environmental data affect the prey (including Anchovy [Engraulus japonicus], Pacific Sand Lance [Ammodytus hexapterus], and Capelin [Mallotus vilosus]) that this species relies on, which will also affect seabird fitness.For example, Pacific Sand Lance are less abundant and grow slower under warmer oceanic conditions (Robards et al., 2002), and embryo survival decreases as pH increases (Murray et al., 2019).

Assis
Given that maximum sea temperature and pH peak during the summer, the variables used in this study reflect environment conditions during the breeding season of rhinoceros auklets.

| Partial redundancy analysis
We used redundancy analysis (RDA) to examine genotype-environment associations for rhinoceros auklets.Redundancy analysis combines ordination and regression techniques to analyze the relationship between multiple predictor variables and a single response variable (Legendre & Legendre, 2012).RDA methods are increasingly being used in landscape genomic studies because this analysis provides a flexible and robust method to analyze the relationship between genetic and environmental variation (Capblancq & Forester, 2021;Forester et al., 2018).Redundancy analysis also allows for variance partitioning to examine the contributions of each variable on the response variable.Another advantage of RDA is that partial RDA models allow for covariates to be included to account for any variables that may influence the relationship between genetic and environmental variation (Forester et al., 2018).For our RDA and partial RDA models, we used individual-based allele count data (genotypes were scored as 0/1/2) for each locus as our full SNP (19,213 SNPs) genetic dataset.
For all missing data, we used an imputation approach where we converted missing genotypes to the most common genotype at each locus across all individuals following the approach of Cayuela et al. (2022).
Our environmental data consisted of the three variables described above: maximum current velocity, pH, and maximum sea surface temperature.Finally, we included population genetic structure as a covariate in our model to account for the signal of population genetic structure on genotype-environment associations using the first two principal coordinates from our population genetic structure analysis with the 9349 SNP dataset for this variable (Forester et al., 2018).We analyzed alternative hypotheses alongside environment.Given that spatial data can have a strong effect on genetic patterns, we analyzed the effect of isolation-by-resistance between the western and eastern Pacific Ocean populations and isolation-by-distance within the eastern and western Pacific Ocean populations.We analyzed the effect of isolation-by-resistance between the eastern and western populations, as opposed to isolation-by-distance, because these two areas are disjunct and separated by deep ocean waters (Hipfner et al., 2020).For this analysis, we used least-cost corridor resistance surface distance data presented by Prill (2019).We converted these data to continuous variables using a principal component analysis in Past 4.11 and used the first and second principal components as predictor variables in our redundancy models.Within the eastern and western Pacific Ocean, we examined isolation-by-distance because Prill (2019) showed that there are no biogeographic barriers preventing individuals from moving between breeding colonies within either the eastern or western Pacific Ocean basins.Therefore, we included latitude and longitude as predictor variables to quantify the effect of geographic distance on genetic patterns.We identified candidate loci using the loading scores for each SNP from the RDA following the approach outlined in Forester et al. (2018).SNPs with loading scores that were ±3 SD (two-tailed p-value = .003)from the mean score of each of the constrained loading scores were considered outlier SNPs.
We used the first three constrained axes as a greater proportion of the variation was explained by these three axes.
To complement our partial RDA, we used the R package latent factor mixed models (LFMM; Caye et al., 2019;Frichot et al., 2013) to identify putative adaptive loci.These models detect correlations between allele frequencies and environmental variation while estimating random effects (including population structure, population history and isolation-by-distance) to reduce the number of false-positive putative loci under selection detected.For these analyses, we used K = 2 based on the overall pattern of population genetic structure observed with PCA and used the same three environmental variables, maximum current velocity, pH, and maximum sea surface temperature, that we used for our partial RDA models.
Once we ran LFMM, we examined the genomic inflation factor (GIF) to assess how well the models accounted for potential confounding factors in the dataset.We then used a false discovery rate of 0.01 to identify candidate loci below this threshold following the approach of François et al. (2016).

| RE SULTS
We  TA B L E 3 Summary table for the redundancy and partial redundancy models examining the effect of environment, population genetic structure, isolation-by-resistance, and isolation-by-distance on genetic variation for rhinoceros auklets.Our model examining genotype-environment associations in the western Pacific accounted for a relatively small portion of the variance (r 2 adj = 0.3%, p = .006).Redundancy models examining population genetic structure and isolation were also significant and accounted for a relatively low proportion of variance (0.1% and 0.2%, respectively).Genetic variation in the western Pacific appears to be associated with spatial data, given that colonies that are geographically closer appear to cluster together (Figure 4a).with the full dataset) and previous studies (Abbott et al., 2014;Hipfner et al., 2020) suggest that deep ocean waters in the Bering Sea act as a barrier to dispersal between these populations (Hipfner et al., 2020;Prill, 2019).Prill (2019)  Maximum (Abbott et al., 2014).
We detected low levels of genetic differentiation among colonies both within the western and Pacific.Philopatry (Coulson, 2016;Friesen, 2015) and spatial segregation during the non-breeding season (Friesen et al., 2007;Munro & Burg, 2017)  One other potential factor that could reduce genetic differentiation among breeding colonies within the eastern and western Pacific is the possibility of sex-linked loci in our dataset.Since we did not have information on the sex of each bird, we were unable to detect and remove putative sex-linked loci from our dataset.
Adaptive candidate loci were associated with the three environmental variables we examined, but a greater proportion were associated with maximum current velocity and pH.Ocean currents and pH play a key role in the distribution and diversity of prey species in marine systems (Ballance et al., 2006;Bost et al., 2009;Chambers et al., 2011;Munro & Burg, 2017;Young et al., 2015).Therefore, the strong association of pH and maximum current velocity with genetic variation could reflect in the differences in foraging behavior and preferences that exist between eastern and western Pacific Ocean populations (Burger et al., 1993;Cunningham et al., 2018;Senzaki et al., 2014;Takahashi et al., 2001;Wilson & Manuwai, 1986).It is well established that when foraging behavior differs between populations, the formation of ecotypes can occur and these ecotypes may become genetically distinct following reduced gene flow (Abeyrama, 2023;Breistein et al., 2022;de Bruyn et al., 2013;Fruet et al., 2017).
Temperature is associated with energy expenditures for many marine species, including rhinoceros auklets (Shimabukuro et al., 2023), yet we found relatively few adaptive candidate loci associated with maximum sea surface temperatures despite western Pacific populations occurring in noticeably warmer waters than eastern Pacific populations.Temperature has been shown to be an important determinant of genetic variation for other marine species (Antoniou et al., 2021;Bradbury et al., 2010;Limborg et al., 2012;Reusch, 2014) and is thought to act as a key driver of diversification for marine species (Bowen et al., 2016), including seabirds (Torres et al., 2021).It is possible that the relationship between temperature and genetic variation may be more important in exothermic organisms that cannot thermoregulate or behaviorally thermoregulate as opposed to endothermic organism like seabirds.The limited influence of temperature on genetic variation detected in this study may also reflect adaptability by rhinoceros auklets.rhinoceros auklets occupy large home ranges (570,000-2,200,000 km 2 ; Hipfner et al., 2020) and inhabit large ranges during the full annual cycle.One of our initial goals was to evaluate this relationship between genetic and sea surface temperature variation, given that sea surface temperatures continue to rise.It is possible that our analysis of this relationship was not sensitive enough (i.e., we did not analyze the whole genome) or other temperaturerelated variables may be more important, and therefore, further analyses are required to characterize this relationship.Additionally, our analyses focused exclusively on breeding colony environmental variables.Previous work on Black-browed Albatross (Thalassarche melanophris), a highly vagile seabird, has suggested that environmental variability on wintering grounds exerts greater effects on fitness than environmental variability during the breeding season (Jenouvrier et al., 2018).Based on the results of this previous study, it is possible that breeding colony environmental variation is a poor predictor of genetic variation for rhinoceros auklets.Further examining the environmental variation that rhinoceros auklets encounter on their breeding grounds does not reflect the true scale of environmental variation that individuals encounter because individuals may potentially spend less time on their breeding grounds than on their wintering grounds.Therefore, selection pressures for adaption outside of the breeding season may exert a strong influence on genetic structure and fitness for rhinoceros auklets as suggested by Jenouvrier et al. (2018).
This study quantifies the relationship between environmental variation and genetic variation for rhinoceros auklets.Our analyses suggest that environmental variation explained a low portion of genetic variation between western and eastern Pacific populations.One possible reason for a small influence of environmental variation on genetic variation is the sequencing method used.
Although RADseq offers a robust method to examine genetic variation, the power of this sequencing method is noticeably less than Figure 3a) when population genetic structure was accounted for.This result indicates that population genetic structure alone does

F
I G U R E 3 (a) Partial redundancy analysis plot for RDA axes 1 and 2 using 19,213 SNPs incorporating the effect of population genetic structure on adaptive genetic variation.Circles (blue) designate individuals sampled from western Pacific breeding colonies, while triangles (light yellow) designate individuals sampled from eastern Pacific breeding colonies.Arrows represent the direction and magnitude of effect by the three environmental variables (maximum current velocity, pH, and maximum sea surface temperature) examined, and loci are shown in gray.(b) Shows the distribution of loci in the ordination space.Candidate loci associated with specific environmental variables are designated by color: purple = maximum current velocity, yellow = pH, and green = maximum sea surface temperature.Venn diagram compares the number of candidate loci detected in the partial RDA model and LFMM.(c) Pie charts show the number of candidate loci that were associated with the three environment variables for the partial RDA model and LFMM.
genetic variation and that genetic-environment associations are present.The three RDA axes explained 67.3%, 17.3%, and 15.4% of the observed variation.Overall, the partial RDA model identified 262 candidate loci that are associated with environmental variation between western and eastern Pacific populations.A greater proportion of loci were associated with pH and maximum current velocity, while fewer loci were associated with maximum sea surface temperature (Figure3b,c).By comparison, the redundancy models examining the influence of population genetic structure (r 2 adj = 4.0%, p = .001)and isolation-by-resistance (r 2 adj = 7.3%, p = .001)on genetic variation accounted for a greater proportion of variance than the partial model that examined the influence of environment.Latent factor mixed models (LFMM) identified 104 candidate loci (p < .01).Ninety-six loci (Figure3b) were identified as candidate adaptive loci by both LFMM and partial RDA model.LFMM identified eight additional candidate adaptive loci that were not identified by the partial RDA model, and more loci were associated with pH (n = 59) and maximum current velocity (n = 37) than with maximum sea surface temperature (n = 8; Figure3c).
The partial redundancy model identified 165 candidate loci, for which the candidate loci and environment associations were evenly distributed across the three environmental variables (maximum current velocity, N = 54; pH, N = 54; and maximum sea surface temperature, N = 57).Our partial redundancy and redundancy models examining genotype-environment associations in the eastern Pacific were not significant (p = .24);environmental variables accounted for a relatively low portion of the observed variance (r 2 adj = 0.1%) among eastern Pacific breeding colonies.The model examining population genetic structure was significant but accounted for a low percentage of variation (r 2 = 0.1%, p = .007).Isolation-by-distance accounted for a slightly higher proportion of variance (r 2 = 0.2%, p = .11),and similar to genetic patterns in the western Pacific Ocean, genetic patterns appear to vary with geographic distance (Figure 4b).The colonies further north, Chowiet and Middleton Islands, clustered further away from the other breeding colonies, while the southernmost colonies, Año Nuevo and southeast Farallon, also clustered separately from all other colonies.4| DISCUSS IONOur genetic analyses of rhinoceros auklet breeding colonies across their range provide evidence for genotype-environment associations.Environmental variation accounted for a relatively low but significant amount of the genetic variation observed between western Pacific and eastern Pacific populations; genetic variation was more strongly associated with pH and maximum current velocity than with maximum sea surface temperature based on partial RDA and LFMM analyses.When we compared our three models (environment, population structure (with 9349 SNPs), and isolation-by-resistance), the isolation-by-resistance model had the highest adjusted-r 2 value, which indicates that physical barriers limit gene flow across the range and have a strong influence on genetic variation for vagile species like seabirds(Lombal et al., 2020).That all three of the models examined are associated with genetic variation reflects that multiple factors drive patterns of genetic variation and differentiation.Overall, we found little evidence of genetic-environment associations within each ocean basin where population genetic patterns appear to follow a pattern of isolation-by-distance.Many seabird populations show reduced population genetic structure over large areas because there are relatively few barriers to gene flow in marine systems(Cowen et al., 2007;Friesen et al., 2007;   Jahnke & Jonsson, 2022).The limited gene flow based on population genetic structure between the western and eastern Pacific F I G U R E 4 (a) RDA plots for RDA axes 1 and 2 using 19,213 SNPs examining the effect of geographic distance on genetic variation in the western Pacific Ocean.(b) RDA plots for RDA axes 1 and 2 using 19,213 SNPs examining the effect of geographic distance on genetic variation in the eastern Pacific Ocean.Arrows represent the direction and magnitude of the effect by latitude and longitude in the model.rhinoceros auklet populations observed in this study (as reflected by the PCoA, DAPC, and pairwise F ST comparisons with the 9349 SNP dataset and the RDA model examining isolation-by-resistance noted that rhinoceros auklets appear to disperse readily between colonies within both the eastern and western Pacific Ocean and therefore hypothesized that other ecological or behavioral factors may act to limit dispersal and gene flow between eastern and western Pacific Ocean populations.The pattern of the western Pacific populations being distinct from the eastern Pacific populations has been observed in some vagile seabird species (González-Jaramillo & Rocha-Olivares, 2011; Hermanet al., 2022), but not all(Pshenichnikova et al., 2015).Population genetic differentiation between the western and eastern Pacific may also reflect the lasting effects of Pleistocene glaciations, where the ranges of the eastern and western Pacific populations shifted southwards isolating them from each other during the Last Glacial are both thought to drive genetic differentiation in seabirds and likely play a role in local genetic patterns within the western and eastern Pacific for rhinoceros auklets.Like other alcids, rhinoceros auklets are known to exhibit strong breeding philopatry(Hipfner et al., 2020;Morrison et al., 2011), and that genetic patterns appear to follow patterns of isolation-by-distance within both ocean basins, suggests that breeding philopatry may explain genetic differences.Non-breeding season distributions also influence genetic patterns.rhinoceros auklets winter at-sea outside of the breeding season, andHipfner et al. (2020) found that genetic differentiation was associated with the extent of spatial overlap during the non-breeding season; populations that exhibited lower spatial overlap were more genetically differentiated.
ST values were considered significantly different from zero at p < .011indicated in bold, following sequential Bonferroni corrections.Colony names are coded and full colony names can be found in Table1.Western Pacific Ocean breeding colonies are highlighted in blue, while eastern Pacific Ocean breeding colonies are highlighted in light yellow.Note: Pseudo-F-score (F), p-value (p), degrees freedom (df), inertia, percent variation explained by each variable (r 2 ), and adjusted percent variation explained by each variable (r 2 adj ) are reported for each model.All bolded values are models considered significant at p < .05.
-Note: Values below diagonal represent pairwise F ST values made among populations with the 9349 SNP dataset.Values above diagonal represent p-values; pairwise F